/*****************************************************************************
 *  ENTROPY - emerging network to reduce orwellian potency yield
 *
 *  Copyright (C) 2002 Juergen Buchmueller <pullmoll@stop1984.com>
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program; if not, write to the Free Software Foundation,
 *  Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA
 *
 *	$Id: fec.c,v 1.2 2005/07/12 23:12:29 pullmoll Exp $
 *****************************************************************************/



/*****************************************************************************
 * This code is based on Luigi Rizzo's sample implementation of an erasure
 * codec. It is based on Vandermode matrices computed over GF(p ^ n).
 * Luigi Rizzo's page is here: http://info.iet.unipi.it/~luigi/fec.html
 * My implementation uses pre-computed tables, generated using Luigi's code.
 * It is specialized and simplified for the MM = 8, n = 8 and idx 0 to 3 case,
 * whereas Luigi's code allows to choose between different values for MM.
 *****************************************************************************/

#include "fec.h"

#define	DEBUG_TABLES	0


#ifdef	TESTFEC
#define	xmalloc(size) malloc(size)
#define	xcalloc(num,size) calloc(num,size)
#define	xfree(ptr) if (NULL != ptr) { free(ptr); ptr = NULL; }
#else
#include "memalloc.h"
#endif

#include "logger.h"

#define MM	8
#define NN	((1<<MM)-1)
#define A0	(NN)	/* log(0) isn't defined, use a special value */
#define	MODNN	(16*(NN+1))

#define N	8

static uint8_t *gfexp = NULL;
static uint8_t *gflog = NULL;
static uint8_t *inverse = NULL;
static uint8_t *gfmult = NULL;

static uint64_t *bitmask = NULL;

#define	GFMULT(i) &gfmult[(i)*(NN+1)]

#undef	MSK
#if	!WORDS_BIGENDIAN
#define	MXOR	0
#else
#define	MXOR	56
#endif

/*
 * Compute x % NN, where NN is 2**MM - 1, without a slow divide.
 * Many many calls, about 1/8 of them is doen with a small x < 2 * NN
 */
static inline uint8_t modnn(uint32_t x)
{
	while (x >= NN) {
		x -= NN;
		x = (x >> MM) + (x & NN);
	}
	return x;
}

/*****************************************************************************
 * fec_encode
 * input:
 *  data    - pointer to 8 * data packets
 *  size    - size of each data packet
 *  idx     - which fec to make
 * output:
 *  fec     - pointer to fec packet
 *****************************************************************************/
int fec_encode(uint8_t *data, size_t size, int idx, uint8_t *fec)
{
	size_t i, j;
	const uint8_t *mul;
	FUN("fec_encode");

	if ((N * idx) >= MODNN) {
		LOG(L_ERROR,("idx out of range (%d)\n", idx));
		errno = EINVAL;
		return -1;
	}
	memset(fec, 0, size);
	for (i = 0; i < N; i++) {
		mul = GFMULT(gfexp[modnn(i * idx)]);
		for (j = 0; j < size; j++)
			fec[j] ^= mul[*data++];
	}
	return 0;
}

/*****************************************************************************
 * decode
 * input:
 *  src     - vector of available data blocks
 *  size    - size of each data block
 *  idx     - indices of available data blocks
 * output:
 *  dst     - reconstructed data, 8 * size bytes
 * returns:
 *  0       - success
 *  -1      - failure to build b matrix or diagonal element is 0
 *****************************************************************************/
#define SWAP(a,b,type) {type tmp; tmp=a; a=b; b=tmp;}
int fec_decode(uint8_t *dst, uint8_t *src[], size_t size, int idx[])
{
	int again = 0;
	int missing = 0;
	uint8_t b[N][N];
	uint8_t a1[N][N];
	const uint8_t *mul;
	int i, j, r, pow, row, col;
	size_t offs;
	FUN("fec_decode");

	/* sort src pointers in ascending idx[] order */
	do {
		again = 0;
		for (i = 0; i < N; i++) {
			j = idx[i];
			if (j < N && j != i) {
				SWAP(src[i], src[j], uint8_t*);
				SWAP(idx[i], idx[j], int);
				again = 1;
				break;
			}
		}
	} while (0 != again);

	/* build b matrix */
	memset(b, 0, sizeof(b));
	memset(a1, 0, sizeof(a1));
	for (i = 0; i < N; i++) {
		b[i][i] = 1;
		if (idx[i] < N) {
			if (i != idx[i]) {
				LOG(L_ERROR,("ouch, %d should not be at %d\n",
					idx[i], i));
				errno = EINVAL;
				return -1;
			}
			a1[i][i] = 1;
		} else {
			pow = idx[i] - N;
			for (j = 0; j < N; j++)
				a1[i][j] = gfexp[modnn(j * pow)];
		}
	}

	/* invert matrix, using a crude method... */
	for (row = 0; row < N; row++) {
		if (0 == a1[row][row]) {
			/* pivot element */
			for (r = row+1; r < N; r++) {
				if (0 != a1[row][r]) {
					SWAP(src[row], src[r], uint8_t *);
					SWAP(idx[row], idx[r], int);
					for (i = 0; i < N; i++) {
						 SWAP(a1[i][row], a1[i][r], uint8_t);
						 SWAP(b[i][row], b[i][r], uint8_t);
					}
					break;
				}
			}
			if (r == N) {
				LOG(L_ERROR,("ouch, diagonal element %d = 0\n",
					 row));
				errno = EINVAL;
				return -1;
			}
		}

		if (1 != (pow = inverse[a1[row][row]])) {
			mul = GFMULT(pow);
			for (col = 0; col < N; col++) {
				b[row][col] = mul[b[row][col]];
				a1[row][col] = mul[a1[row][col]];
			}
		}

		for (r = 0; r < N; r++) {
			if (r == row || 0 == (pow = a1[r][row]))
				continue;
			mul = GFMULT(pow);
			if (row == idx[row]) {
				/* source, only a1[row][row] != 0 */
				b[r][row] ^= mul[b[row][row]];
				a1[r][row] ^= mul[a1[row][row]];
				continue;
			}
			for (col = 0; col < N; col++) {
				b[r][col] ^= mul[b[row][col]];
				a1[r][col] ^= mul[a1[row][col]];
			}
		}
	}

	/* now comes the actual decoding of missing blocks */
	for (row = 0; row < N; row++) {

		if (row == idx[row]) {
			memcpy(&dst[row * size], src[row], size);
			continue;
		}

		missing++;
		memset(&dst[row * size], 0, size);
		for (col = 0; col < N; col++) {
			uint8_t *d = &dst[row * size];
			uint8_t *s = src[col];
			mul = GFMULT(b[row][col]);

			if (0 == (size % 8)) {
				uint32_t *dd = (uint32_t *)d;
				for (offs = 0; offs < size; offs += 8) {
#if	!WORDS_BIGENDIAN
					register uint32_t l = *(uint32_t *)&s[offs+0];
					register uint32_t h = *(uint32_t *)&s[offs+4];
					*dd++ ^=
						((uint32_t)mul[(uint8_t)(l >>  0)] <<  0) |
						((uint32_t)mul[(uint8_t)(l >>  8)] <<  8) |
						((uint32_t)mul[(uint8_t)(l >> 16)] << 16) |
						((uint32_t)mul[(uint8_t)(l >> 24)] << 24);
					*dd++ ^=
						((uint32_t)mul[(uint8_t)(h >>  0)] <<  0) |
						((uint32_t)mul[(uint8_t)(h >>  8)] <<  8) |
						((uint32_t)mul[(uint8_t)(h >> 16)] << 16) |
						((uint32_t)mul[(uint8_t)(h >> 24)] << 24);
#else
					register uint32_t h = *(uint32_t *)&s[offs];
					register uint32_t l = *(uint32_t *)&s[offs+4];
					*dd++ ^=
						((uint32_t)mul[(uint8_t)(h >> 24)] << 24) |
						((uint32_t)mul[(uint8_t)(h >> 16)] << 16) |
						((uint32_t)mul[(uint8_t)(h >>  8)] <<  8) |
						((uint32_t)mul[(uint8_t)(h >>  0)] <<  0);
					*dd++ ^=
						((uint32_t)mul[(uint8_t)(l >> 24)] << 24) |
						((uint32_t)mul[(uint8_t)(l >> 16)] << 16) |
						((uint32_t)mul[(uint8_t)(l >>  8)] <<  8) |
						((uint32_t)mul[(uint8_t)(l >>  0)] <<  0);
#endif
				}
			} else for (offs = 0; offs < size; offs++)
				*d++ ^= mul[*s++];
		}
	}

	return 0;
}

/*****************************************************************************
 * fec84_encode
 * input:
 *  in      - pointer to size data bytes
 *  size    - number of data bytes
 *	flags	- 12 flags telling which bits are already there (1 bits)
 * output:
 *  out     - pointer to buffer receiving 12 chunks
 *****************************************************************************/
int fec84_encode(uint8_t *out, uint8_t *in, size_t size, int flags)
{
	register size_t offs, hsize = (size + 7) / 8;

	/* return zero if all twelve bits are already there */
	if (0xfff == flags) {
		return 0;
	}

	/* split into bits only if at least one column is missing */
	if (0xff != (flags & 0xff)) {
		for (offs = 0; offs < hsize; offs++, in += 8) {
			register uint64_t w = *(uint64_t *)in;
			register uint64_t bswap =
				bitmask[0*256+(uint8_t)(w >> ( 0^MXOR))] |
				bitmask[1*256+(uint8_t)(w >> ( 8^MXOR))] |
				bitmask[2*256+(uint8_t)(w >> (16^MXOR))] |
				bitmask[3*256+(uint8_t)(w >> (24^MXOR))] |
				bitmask[4*256+(uint8_t)(w >> (32^MXOR))] |
				bitmask[5*256+(uint8_t)(w >> (40^MXOR))] |
				bitmask[6*256+(uint8_t)(w >> (48^MXOR))] |
				bitmask[7*256+(uint8_t)(w >> (56^MXOR))];
			if (0 == (flags & 0x001))
				out[0*hsize+offs] = (uint8_t)(bswap >> ( 0^MXOR));
			if (0 == (flags & 0x002))
				out[1*hsize+offs] = (uint8_t)(bswap >> ( 8^MXOR));
			if (0 == (flags & 0x004))
				out[2*hsize+offs] = (uint8_t)(bswap >> (16^MXOR));
			if (0 == (flags & 0x008))
				out[3*hsize+offs] = (uint8_t)(bswap >> (24^MXOR));
			if (0 == (flags & 0x010))
				out[4*hsize+offs] = (uint8_t)(bswap >> (32^MXOR));
			if (0 == (flags & 0x020))
				out[5*hsize+offs] = (uint8_t)(bswap >> (40^MXOR));
			if (0 == (flags & 0x040))
				out[6*hsize+offs] = (uint8_t)(bswap >> (48^MXOR));
			if (0 == (flags & 0x080))
				out[7*hsize+offs] = (uint8_t)(bswap >> (56^MXOR));
		}
	}

	if (0 == (flags & 0x100))
		fec_encode(out, hsize, 0, &out[ 8*hsize]);
	if (0 == (flags & 0x200))
		fec_encode(out, hsize, 1, &out[ 9*hsize]);
	if (0 == (flags & 0x400))
		fec_encode(out, hsize, 2, &out[10*hsize]);
	if (0 == (flags & 0x800))
		fec_encode(out, hsize, 3, &out[11*hsize]);

	return 0;
}

/*****************************************************************************
 * fec84_decode
 * input:
 *  in      - pointer to data blocks
 *  size    - size of resulting data
 *  flags   - flags indicating available blocks
 * output:
 *  out     - pointer to resulting data
 *****************************************************************************/
int fec84_decode(uint8_t *out, uint8_t *in, size_t size, int flags)
{
	uint64_t *dq = (uint64_t *)out;
	register size_t offs, hsize = (size + 7) / 8;
	uint8_t *src[N];
	int i, j, idx[N];
	int rc;
	FUN("fec84_decode");

	/* construct index array from available blocks */
	for (i = 0, j = 0; i < 12; i++) {
		if (0 == (flags & (1 << i)))
			continue;
		idx[j] = i;
		src[j] = &in[i * hsize];
		if (++j >= N)
			break;
	}

	/* we cannot reconstruct data with less than 8 blocks */
	if (j < N) {
		errno = EINVAL;
		rc = -1;
		return rc;
	}

	if (0 != (rc = fec_decode(in, src, hsize, idx))) {
		LOG(L_ERROR,("fec_decode() call failed\n"));
		return rc;
	}

	for (offs = 0; offs < hsize; offs++) {
		register uint64_t bswap =
			bitmask[0*256+in[0*hsize+offs]] |
			bitmask[1*256+in[1*hsize+offs]] |
			bitmask[2*256+in[2*hsize+offs]] |
			bitmask[3*256+in[3*hsize+offs]] |
			bitmask[4*256+in[4*hsize+offs]] |
			bitmask[5*256+in[5*hsize+offs]] |
			bitmask[6*256+in[6*hsize+offs]] |
			bitmask[7*256+in[7*hsize+offs]];
		*dq++ = bswap;
	}

	return 0;
}

#undef	SIZE
#define SIZE 4096

int fec(void)
{
	/* primitive polynomial of GF(2) for MM = 8 */
	/* 1 + x^2 + x^3 + x^4 + x^8 */
	static const char Pp[] = "101110001";
	int flags;
	uint8_t *raw = NULL, *enc = NULL, *ecp = NULL, *cmp = NULL;
	uint32_t mask;
	uint32_t a, b, c, d, i, j;
	uint64_t m;
	time_t t0, t;
	FUN("fec");

#if	DEBUG_TABLES
	printf("bitmask[%d*%d]\n", 8, 256);
#endif
	bitmask = (uint64_t *)xcalloc(8 * 256, sizeof(uint64_t));
	for (i = 0; i < 8; i++) {
#if	DEBUG_TABLES
		printf("# shift %d\n", i);
#endif
		for (j = 0; j < 256; j++) {
			m = 0;
#if	!WORDS_BIGENDIAN
			if (j & 0x01) m |= (0x0000000000000001llu << i);
			if (j & 0x02) m |= (0x0000000000000100llu << i);
			if (j & 0x04) m |= (0x0000000000010000llu << i);
			if (j & 0x08) m |= (0x0000000001000000llu << i);
			if (j & 0x10) m |= (0x0000000100000000llu << i);
			if (j & 0x20) m |= (0x0000010000000000llu << i);
			if (j & 0x40) m |= (0x0001000000000000llu << i);
			if (j & 0x80) m |= (0x0100000000000000llu << i);
#else
			if (j & 0x01) m |= (0x0100000000000000llu << i);
			if (j & 0x02) m |= (0x0001000000000000llu << i);
			if (j & 0x04) m |= (0x0000010000000000llu << i);
			if (j & 0x08) m |= (0x0000000100000000llu << i);
			if (j & 0x10) m |= (0x0000000001000000llu << i);
			if (j & 0x20) m |= (0x0000000000010000llu << i);
			if (j & 0x40) m |= (0x0000000000000100llu << i);
			if (j & 0x80) m |= (0x0000000000000001llu << i);
#endif
			bitmask[i * 256 + j] = m;
#if	DEBUG_TABLES
			printf("%016llx,", bitmask[i * 256 + j]);
			if (3 == (j % 4))
				printf("\n");
#endif
		}
	}

	/* initialize index -> polynomial and polynomial -> index lookup tables */
	gfexp = (uint8_t *)xcalloc((NN+1), sizeof(uint8_t));
	gflog = (uint8_t *)xcalloc((NN+1), sizeof(uint8_t));
	for (i = 0, mask = 1; i < MM; i++) {
		gfexp[i] = mask;
		gflog[gfexp[i]] = i;
		/* if Pp[i] == 1 then term @^i occurs in poly repr of @^MM */
		if (Pp[i] == '1')
			gfexp[MM] ^= mask;
		mask <<= 1;
	}
	gflog[gfexp[MM]] = MM;

	/*
	 * Now we have the poly representation of @^MM.
	 * The poly representation of @^(i+1) is given by poly repr. of
	 * @^i shifted left by one bit and accounting for any @^MM term
	 * that may occur when poly repr. of @^i is shifted.
	 */
	mask >>= 1;
	for (i = MM + 1; i < NN; i++) {
		if (gfexp[i-1] >= mask)
			gfexp[i] = gfexp[MM] ^ ((gfexp[i-1] ^ mask) << 1);
		else
			gfexp[i] = gfexp[i-1] << 1;
		gflog[gfexp[i]] = i;
	}
	gflog[0] = A0;
	gfexp[NN] = 0;

#if	DEBUG_TABLES
	printf("gfexp[%d]\n", NN+1);
	for (i = 0; i < NN+1; i++) {
		printf("%3d,", gfexp[i]);
		if (15 == (i % 16))
			printf("\n");
	}
	printf("gflog[%d]\n", NN+1);
	for (i = 0; i < NN+1; i++) {
		printf("%3d,", gflog[i]);
		if (15 == (i % 16))
			printf("\n");
	}
#endif
	/* inverse of a number inv[@^i] = @^[NN-i-1] */
	inverse = (uint8_t *)xcalloc(NN + 1, sizeof(uint8_t));
	inverse[0] = NN;	/* invalid! */
	inverse[1] = 1;
	for (i = 2; i < NN+1; i++)
		inverse[i] = gfexp[NN - gflog[i]];

#if	DEBUG_TABLES
	printf("inverse[%d]\n", NN+1);
	for (i = 0; i < NN+1; i++) {
		printf("%3d,", inverse[i]);
		if (15 == (i % 16))
			printf("\n");
	}
#endif

	/* initialize a multiplication lookupt table for two numbers modulo NN */
	gfmult = (uint8_t *)xcalloc((NN+1) * (NN+1), sizeof(uint8_t));
	for (i = 0; i < NN+1; i++)
		for (j = 0; j < NN+1; j++)
			gfmult[i*(NN+1)+j] = gfexp[modnn(gflog[i] + gflog[j])];
	for (j = 0; j < NN+1; j++)
		gfmult[0*(NN+1)+j] = gfmult[j*(NN+1)+0] = 0;

#if	DEBUG_TABLES
	printf("gfmult[%d*%d]\n", NN+1, NN+1);
	for (i = 0; i < NN+1; i++) {
		printf("# row %3d\n", i);
		for (j = 0; j < NN+1; j++) {
			printf("%3d,", gfmult[i*(NN+1)+j]);
			if (15 == (j % 16))
				printf("\n");
		}
		printf("\n");
	}
	return -1;
#endif

#if	TESTFEC == 0
	info("tables ");
#endif

	/* Test run of the FEC codec */
	raw = (uint8_t *)xcalloc(8 * SIZE, sizeof(uint8_t));
	enc = (uint8_t *)xcalloc(12 * SIZE, sizeof(uint8_t));
	ecp = (uint8_t *)xcalloc(12 * SIZE, sizeof(uint8_t));
	cmp = (uint8_t *)xcalloc(8 * SIZE, sizeof(uint8_t));

	/* Fill the buffer 'raw' with 8*SIZE bytes of pseudo random data */
	srand(time(NULL));
	for (i = 0; i < 8 * SIZE; i++)
		raw[i] = (rand() >> 16) % 256;
	if (0 != fec84_encode(enc, raw, 8 * SIZE, 0)) {
		LOG(L_ERROR,("fec84_encode() call failed\n"));
		die(1, "fec84_encode() call failed\n");
	}
	memcpy(ecp, enc, 12 * SIZE);
	t0 = time(NULL);
	/* check all combinations of four erasures */
	for (a = 0; a < 9; a++) {
		for (b = a+1; b < 10; b++) {
			for (c = b+1; c < 11; c++) {
				for (d = c+1; d < 12; d++) {
					t = time(NULL);
#if	TESTFEC == 0
					if (t != t0) {
						info("%x:%x:%x:%x\b\b\b\b\b\b\b", a, b, c, 11);
						t0 = t;
					}
#endif
					flags = 0;
					for (i = 0; i < 12; i++) {
						if (i == a || i == b || i == c || i == d) {
							memset(&enc[i * SIZE], 0, SIZE);
						} else {
							memcpy(&enc[i * SIZE], &ecp[i * SIZE], SIZE);
							flags |= 1 << i;
						}
					}
					/* now decode this erasure combination */
					if (0 != fec84_decode(cmp, enc, 8*SIZE, flags)) {
						LOG(L_ERROR,("fec84_decode(%x:%x:%x:%x) call failed\n",
							a, b, c, d));
						die(1, "fec84_decode(%x:%x:%x:%x) call failed\n",
							a, b, c, d);
					}
					/* data must match now */
					if (0 != memcmp(raw, cmp, 8*SIZE)) {
						LOG(L_ERROR,("memcmp(%x:%x:%x:%x) call failed\n",
							a, b, c, d));
						die(1, "memcmp(%x:%x:%x:%x) call failed\n",
							a, b, c, d);
					}
				}
			}
		}
	}
	xfree(raw);
	xfree(enc);
	xfree(ecp);
	xfree(cmp);

	info("checked ");
	return 0;
}


#ifdef	TESTFEC
char progname[] = "fec";
configuration_t *g_conf = NULL;

void hexdump(uint8_t *data, size_t size)
{
	size_t i, j;
	for (i = 0; i < size; i++) {
		if (0 == (i % 16)) {
			printf("%04x:", i);
		}
		printf(" %02x", data[i]);
		if (15 == (i % 16)) {
			printf(" - ");
			for (j = i - (i % 16); j <= i; j++) {
				int c = data[j];
				printf("%c", c < 32 || c > 126 ? '.' : c);
			}
			printf("\n");
		}
	}
	if (0 != (i % 16)) {
		for (j = i % 16; j < 16; j++) {
			printf("   ");
		}
		printf(" - ");
		for (j = i - (i % 16); j < i; j++) {
			int c = data[j];
			printf("%c", c < 32 || c > 126 ? '.' : c);
		}
		printf("\n");
	}
}

int info(const char *fmt, ...)
{
	int len;

	va_list ap;
	va_start(ap, fmt);
	len = vfprintf(stdout, fmt, ap);
	va_end(ap);
	return len;
}

int die(int rc, const char *fmt, ...)
{
	va_list ap;
	va_start(ap, fmt);
	vfprintf(stderr, fmt, ap);
	va_end(ap);
	exit(rc);

	/* not reached */
	return 0;
}

#undef	SIZE
#define SIZE 16384
int main(int argc, char *argv[])
{
	FILE *fp = NULL;
	uint8_t raw[8*SIZE];
	uint8_t enc[12*SIZE];
	uint8_t ecp[12*SIZE];	/* encode-copy  */
	uint8_t dec[8*SIZE];
	int flags, k, l, m, n, t ;
	size_t i;
	struct timeval tstart, tstop;
	uint64_t bytes, time, rate;
	int rc, cnt;

	(void)argc;
	(void)argv;

	g_conf = xcalloc(1, sizeof(configuration_t));
	strcpy(g_conf->logfile, "fec.log");
	for (i = 0; i < LOGSECTIONS; i++)
		g_conf->loglevel[i] = L_DEBUG;

	if (0 != fec()) {
		fprintf(stderr, "fec() call failed\n");
		exit(1);
	}

	fp = fopen("/dev/urandom", "r");
	if (NULL == fp) {
		perror("/dev/urandom");
		exit(1);
	}
	printf("Tests pm-fec84:\n");
	system("date");
	for (t = 0; t < 20; t++) {		/* total tests  */
		if (8*SIZE != fread(raw, 1, 8*SIZE, fp)) {
			perror("fread(\"/dev/urandom\")");
			exit(2);
		}
		rc = fec84_encode(enc, raw, 8*SIZE);
		if (0 != rc) {
			printf("****** encode error %d:\n", rc);
			exit(3);
		}
		memcpy(ecp, enc, 12*SIZE);

		cnt = 0;
		gettimeofday(&tstart, NULL);
		for (k = 0; k < 9; k++) {
			for (l = k+1; l < 10; l++) {
				for (m = l+1; m < 11; m++) {
					for (n = m+1; n < 12; n++) {
						memcpy(enc, ecp, 12*SIZE);

						flags = 0xfff;
						flags &= ~(1 << k);
						flags &= ~(1 << l);
						flags &= ~(1 << m);
						flags &= ~(1 << n);

						cnt++;
						for (i = 0; i < 12; i++) 
							if (0 == (flags & (1 << i))) 
								memset(&enc[i*SIZE], 0, SIZE);
						rc = fec84_decode(dec, enc, 8*SIZE, flags);
						if (0 != rc) {
							printf( "**** decode error %d:\n", rc);
							hexdump(dec, 8*SIZE);
							exit(4);
						}
						if (0 != memcmp(dec, raw, 8*SIZE)) {
							printf( "**** bad product, src is:\n" );
							hexdump(raw, 8*SIZE);
							rc = 42;
						}
						if (0 != rc) {
							printf( "decoded\n" );
							hexdump(dec, 8*SIZE);
							printf("\nDecode-ERROR %d, block %d, erasure [%d,%d,%d,%d]", 
								rc, t, k, l, m, n );
							exit(5);
						}
					}
				}
			}
		}
		gettimeofday(&tstop, NULL);
		bytes = cnt * 8 * SIZE;
		time = (uint64_t)1000000 * (tstop.tv_sec - tstart.tv_sec) +
			tstop.tv_usec - tstart.tv_usec;
		rate = (uint64_t)1000000 * bytes / time;
		printf( "test # %d Ok! cnt = %d (%lld bytes), time %lldus, %lld MB/s\n",
			t, cnt, bytes, time, rate/1024/1024);
	}
	system("date");
	return rc;
}
#endif	/* TESTFEC */
